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The problem of electromagnetic radiation traveling in a general 
homogeneous, but anisotropic and gyrotropic, medium has been 
solved. The plane wave representation is used to convert Maxwell's 
equations into a general eigenproblem, which allows for tensor 
permeability and permittivity and for electric and magnetic gyro- 
tropy. This formulation can be applied directly to computer evalua- 
tion of wave velocities and polarizations for a given wave-normal 
direction. The reflection and refraction directions and amplitudes at 
an interface between two different, general anisotropic and gyrotropic 
media are solved. It is outlined how, with the aid of a digital 
computer, ray tracing through a general succession of anisotropic 
and gyrotropic, but locally homogeneous, media can be carried out. 

I. INTRODUCTION 

Traditionally, the theory of light propagation in transparent, aniso- 
tropic media has assumed that the permeability of the medium was 
the same as that of the vacuum, or at least that it was isotropic. 
Conversely, the theory of microwave propagation in ferrites and gar- 
nets has assumed that the permittivity was isotropic in spite of the 
permeability being anisotropic and possibly, gyrotropic. (References 1 
and 2 are examples of the optics tradition and the microwave tradition, 
respectively.) 

Both viewpoints are incomplete, particularly, when applied to in- 
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frared optics where anisotropy effects in both permeability and per- 
mittivity could be important. As a case in point, Yttrium Iron Garnet 
(yig) is optically transparent at infrared wavelengths beyond 1 jum 3 
and at microwave frequencies; it could be used at wavelengths where 
the effects of permeability and permittivity are comparable. 

We develop the theory of electromagnetic wave propagation in a 
homogeneous medium where both the permeability and permittivity 
are anisotropic. Furthermore, by assuming that the permittivity and 
permeability tensors are Hermitian, not merely real and symmetric, 
we incorporate the more general effects of gyrotropy. 

II. BASIC ASSUMPTIONS 

We assume the general validity of Maxwell's equations without 
"source terms" (no charges or currents): 

V X H = dD/dt (1) 

V X E = -BB/dt (2) 

V-D = (3) 

V-B = 0, (4) 

where E, D, B, H are complex vector fields which are functions of 
space and time (r, t). 

We shall restrict our considerations to plane-wave solutions where 
the above four vector fields can be cast in the form 

F(r, t) = Foexp[i(k-r - «*)], (5) 

where Fo is a complex vector constant and u is a real constant. The 
vector, k, is usually considered to be real, but may be complex even in 
lossless media (for "evanescent waves"). 

To complete the specifications, we need to define material-depend- 
ent relationships between the four vector fields of Maxwell's equations. 
We do so by means of the electromagnetic energy density: 

I/=%E*.D + %H*.B. (6) 

For a plane wave, we can substitute (5) into (6) to show that U is 
independent of time and position (assuming k is real): 

U m ^Eo* -Do + WHS -Bo. (7) 

By assuming that Uis a purely quadratic function of the components 
of E* , we obtain 

Do - eE . (8) 

The Hermiticity of ?, the permittivity or dielectric tensor, follows from 
the reality of U. 
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Similarly, by assuming that U is a purely quadratic function of the 
components of Ho, we conclude that 

Bo = pHo, (9) 

where the permeability tensor,/!, is Hermitian. 

Our assumptions are consistent with magneto-optic effects where, 
conventionally, a Polder tensor 4 is defined, e.g., 

/ H ifia 0\ 
iL=\-iHa I* . (10) 

V ii' J 

Using the plane-wave representation (5) in Maxwell's first pair of 
equations, (1) and (2), we obtain 

k X H = -coD. (11) 

k X E = wB. (12) 

Similar substitutions into Maxwell's last pair of equations, (3) and 
(4), 

k-D = (13) 

k-B = 0, (14) 

which, in fact, follow from (11) and (12) when w#0. 

III. THE WAVE VELOCITY 

We rewrite the plane-wave representation (5) in a suggestive form 
F = F exp[-(Imk).r]-exp[i|Re k|(n-r - v w t)], (15) 

where 

n = Rek/|Rek| (16) 

is a unit vector that is normal to the wave fronts, and the wave velocity 

Vo, = w/|Rek|. (17) 

Thus, Vu, represents the rate at which the wave fronts appear to be 
advancing in the n direction. 

We usually assume that Im k = 0. However, for Im k ^ we have 
an "evanescent wave"; a wave whose intensity diminishes in the Im k 
direction. It will be shown later, under the discussion of Poynting's 
vector, that Im k is perpendicular to Poynting's vector in nondissipa- 
tive media. 

Let us assume now that Im k = and that, consequently, n is 
parallel to k. We can restate (25) and (26) as follows: 

n X H = -v„,D = -v u ,f E (18) 
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nxE = v„B = v^jBH. (19) 

If we go to a six-dimensional representation, we may combine (18) 
and (19) into a single-matrix equation: 




(20) 



For convenience, we shall compress (20) by defining a more compact 
notation 

or 

av = Wyfiv. (22) 

Note that (22) is clearly a general eigenproblem. The tensor, % is 
determined by the permeability and permittivity of the medium. If k 
is real, then a is determined completely by n. 

Both 5 and j& are usually Hermitian. Tensor, a, is symmetric. It is 
real if k was real. We can restate (6) as 

2U=v*.$v. (23) 

Away from dramatic resonances, we may assume that U is positive 
definite as well as real. If so, j& is not only Hermitian but positive 
definite. 

When ]& is both Hermitian and positive definite, we can find its 
inverse and its square root (all of which are Hermitian and positive 
definite as well). Thus, we are able to convert (8) to an ordinary 
eigenproblem of a Hermitian matrix: 

$-W0-i/*)$V* v ) = v w 1/2 v). (24) 

The Hermiticity of ($~ 1/2 a$" 1/2 ) guarantees that there are six eigen- 
solutions, each with a real v„,. 
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We can use the Hermiticity of 3 and /? to prove orthogonality: 

p*-av2 ss v*i>»* >flv2 = Vuavf -flpz. (25) 

Clearly, when \t,\ does not equal v^, the corresponding eigenvectors 
satisfy the orthogonality relations: 

pf.aP2 = 0= pf-$p 2 . (26) 

Additional orthogonality relationships can be proved if v£,i t* vJk 
(71 = 1,2,3...): 

pf>#($- 1 a) n p 2 = 0. (27) 

Two eigensolutions can be provided automatically: 

" 5 = (oj *** v « = \l 

v„, 5 = and v^ = 0. (28) 

Thus, we have to find four more eigensolutions to this six-dimensional 
eigenproblem. 

The reality of n and the special form of $ enable us to simplify the 
problem still further. Suppose we have a solution to (21). It follows 
then that 

ft 1)(4)-- G >)(-«)• 

Thus, for each eigensolution with wave velocity v w , there is an analo- 
gous solution that travels in the opposite direction. 

In summary, of the six solutions to the eigenproblem (20) for a 
specified, real n, two are trivial and have zero wave velocity, two are 
positive and correspond to two polarized solutions propagating in the 
+n direction, and the last two are negative, and, consequently, travel 
in the — n direction. 

It is worthwhile noting that we have included nonreciprocal cases, 
since Faraday rotation, the traditional nonreciprocal effect, is mani- 
fested by complex, yet Hermitian, permittivity and/or permeability 
tensors. These effects do not seem to cause the wave velocities in the 
— n direction to differ from those in the +n direction. The nonreciprocal 
effects must be manifest in the differing polarizations between the 
forward and the backward directions. 

Polarization is specified usually by the orientation of the electro- 
magnetic field components normal to the relevant velocity direction 
(here n). Thus, from eqs. (13) and (14) we should concentrate upon D, 
when k is real. In general, the components of D, say, are complex; then 
the polarization should be considered elliptical. The special case of 
D oc D* is called "linear polarization." Since D is confined to the plane 
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perpendicular to k, its description for a specified n can be reduced 
always to two complex vector components; this two-dimensional rep- 
resentation is known as the Jones vector. 5,6 

Nonreciprocity is the inequivalence of wave velocity (or attenuation, 
which we ignore here) for a given polarization for the ±n directions. 
The reason that Faraday rotation is nonreciprocal is that the polari- 
zation for a given magnitude of wave velocity in the ±n directions are 
not the same; the Jones vector representation for the ±n directions 
associate right-handed elliptical polarization to left-hand elliptical 
polarization in the opposite direction. Two polarizations in the ±n 
directions would be considered the same if the corresponding D fields 
obeyed the proportionality 

Di oc D 2 *. (30) 

Thus, nonreciprocity for solutions consistent with (20) could occur 
only when D* is not proportional to D; no nonreciprocity occurs for 
modes that are linearly polarized. It also appears to be a consequence 
of the form of (20) that no polarization-independent nonreciprocity 
can be produced in homogeneous media. However, if the medium is 
bianisotropic, then polarization-independent nonreciprocity is possible 
(see Ref. 7). 

The orthogonality relations (26) and (27) derived earlier, can be 
used to show that the two polarizations associated with the two 
different wave velocities in a given ra-direction are "perpendicular." 
Let us call these solutions, "fast" and "slow." By assuming that co ^ 0, 
we can derive from the orthogonality relation 

ELt-Ddow + HLt-Bdow = 0. (31) 

Since EfaBt, — Hfast is again a solution with yet another v w , we can derive 
from orthogonality 

EfasfDalow — Hfast'Bslow = 0. (32) 

Thus, we can conclude from (31) and (32) that 

ELfDslow = (33) 

HLt-Bsiow = 0. (34) 

We see from (31) also that there is no interference in the energy 
density between the fast and slow solutions for a given direction, n. 
Suppose at some location 

V = aVfast + ftPalow, (35) 

then 

U = y*v>$v = K|a| s riU-iftft H t + fc|ft|VU-jkb». (36) 
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Usually, the literature considers the more restrictive cases, 8 where 
either /t or if is isotropic and conclude 

DLfD 8 , ow = 0(?), 

which is not always true for more general cases such as are considered 
here. 

IV. THE POYNTING VECTOR 
If we consider the time dependence of U, we derive: 

dU ldr*j8 1 lk <&r n ( . d?v\ 
dt 2 dt 2 dt \ dt ) 

By using MaxweU's eqs. (1) and (2), we obtain further: 

dU/dt = Re(E*-V x H - H*V x E) = V-[Re(E* x H)]. (38) 

As we have identified dU/dt with the divergence of a "flux," this flux 
should be identified with Poynting's vector. With this in mind, we 
make two definitions 

G = E*XH (39) 

and 

S = Re G. (40) 

For a plane wave of real k, G is independent of time and space. 
However, it may be complex. The real part of G, namely, S, is 
Poynting's vector and is always real, obviously. 

The vectors G and G* are useful because they are "dual" to k. 
Consider 

G x D = H(E*.D) - E*(HD) 

= H(E* D) + w _1 E*(H-k x H) = t/H. 

Consequently, 

G X D = UK. (41) 

A similar derivation shows 

G*xB = -£/E. (42) 

These two equations are analogous to (12) and (11) and suggest a 
series of steps that are usable for calculating the ray velocity. 

The ray velocity vector, v r , which is parallel to the Poynting vector, 
S, is defined as 

\ r =S/U (43a) 

v r =|v r |. (43b) 
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The Appendix proves that the Poynting vector is indeed the ray 
direction. 
Unlike the theory for k, we will not be able to assume that G is real. 

Nevertheless, we define 

£^G/|ReG| (44) 

and 

s = Re#=S/|S|. (45) 

With these notational conventions, we restate (5) and (6) 

g*xB = -v7 a E = -v" 1 ?^ (46) 

g X D = v^H = v-^B. (47) 

These two relations can be rewritten as an eigencondition: 
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= v r ' 1 ?- 1 %* l (48) 



And a terser notation, analogous to (22), could be 

Tv = v^-'v, (49) 

where V is Hermitian; /§ is positive definite and Hermitian because 

2U=v*-$- 1 v (50) 

in analogy to (23). 

In accordance with the analogous discussion of the solutions (22), 
the six real eigenvalue solutions to (49) will consist of two zero v, \ 
two positive v^ 1 ', and their two negatives. Similar orthogonality rela- 
tions can also be generated. 

The duality 9 would be complete if g were equal to s. If k and $ were 
real, we could expect that G = S. However, we have demonstrated by 
computer that a real k does not always lead to a real G. 

Thus, we must conclude that to solve for v r starting with a given s, 
we must find the appropriate values of Im g that produce a real k in 
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a test such as 

£/k = w(D*XB). (51) 

We shall show now that if k is complex, then its imaginary part is 
perpendicular to S. From (11) and (12), we can generate 

-E*kxH = <oE*D (52) 

H*kxE = wH*B. (53) 

Adding (52) and (53), we obtain 

k-(E* X H + E x H*) = <o(E*.D + H*B). (54) 

Using the definitions of Poynting's vector and energy density, we can 
restate (54) as 

k-S = uU. (55) 

Because to, U, and S are always assumed to be real, we can conclude 
that 

S-Imk = 0. (56) 

V. REFLECTION AND REFRACTION 

We can find the reflection and refraction of a plane wave at a flat 
interface between two different, homogeneous media of the general 
type we have been considering. 

Figure 1 indicates what we should expect at the interface between 
two media. The ko is the propagation vector of the incident wave. 
Usually, there will be two reflected waves, represented by ki and k 2 , 
and two refracted (transmitted) waves, represented by k3 and k^ 
Details of the four "scattered" waves are determined by boundary 
conditions. 

To determine ki, k 2 , k3, and ki for a particular ko, the boundary 
conditions we need to consider are the equality of the co(k) associated 
with all five k vectors. In addition, because the geometry possesses 
two-dimensional translational symmetry in the plane of the interface, 
we must match the components in the plane of the interface for all the 
k vectors. These constraints are sufficient to determine the four 
"scattered" k vectors for a given k . 

Finding the scattered k vectors amounts to specifying the compo- 
nents of y = k/co that are in the plane of the interface, yy, and solving 
an eigenproblem that provides 

yl 1 ■«/!!*.£. | . (57) 

We define x ± , for notational convenience, as a unit vector perpendicular 
to the interface and directed in the same sense as the incident radiation. 
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Fig. 1— Scattering at an interface between two anisotropic media. 

From (11) and (12), we have 

-yXH = ?E 
YXE = j2H. 

Rearranging, we have 

-yjfcj. X H = eE + Yh X H 

yjLxE-pH- Y||XE. 

Thus, we have a general eigenproblem again: 



(58) 
(59) 

(60) 
(61) 




= a^v = y 



(62) 
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This eigenproblem has a lot of similarity to (20) in that the matrices 
on each side of the equation are Hermitian. Further, we expect two 
extraneous solutions for which the eigenvalues yl 1 = and the eigen- 
vectors 

*-(jn and P6== {°xX (63) 

are analogous to (28). We expect two "reflected" solutions with nega- 
tive 

§j. = Xj.'S. (64) 

The sign of y± is not so relevant because k and S are rarely parallel in 
anisotropic media. 

We expect a maximum of two "refracted" solutions when we use ?2 
and p 2 in (62) and look for solutions with positive Sj_. We shall have 
less than the maximum number of refracted solutions when yj is getting 
large enough so that /L ceases to be positive definite. 

If /L is neither positive definite nor negative definite, a derivation 
similar to (24) does not exist and some of the solutions for yl 1 may be 
complex. If yj. is complex, then the corresponding E and H fields 
satisfy from (26): 

r*.a^-0. (65) 

Therefore, we have an evanescent solution: 

&. - 0. (66) 

Note that the sign of Yj. is not the criterion for acceptable solutions 
to (62), rather the sign of Sj. is the criterion. 

VI. SCATTERING AMPLITUDES 

The solution of (62) has provided the values of ki, k 2 , fo, and ku for 
a given ko. Let us label the corresponding eigenvectors: v\, V2, Pa, V4, 
and vo. We shall assume no particular normalization for the eigenvec- 
tors, but suppose for a given vo there will be particular amplitudes am 
for each scattered wave. Naturally, the solutions for the a, are depend- 
ent upon the normalizations of the Pi. 

From Maxwell's eqs. (1) and (2), we must satisfy continuity of the 
components of E and H in the plane of the interface. These conditions 
may be restated as 

So-Po = 3j_(— a\V\ — 02^2 + a^Pz + gup-i) • (67) 

Equation (67) can be restated as four simultaneous equations by 
suitable definitions. In the following, assume the index, i, takes on the 
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integer values 1 through 4. 


Let 




bt ■ = v?-cL.vo 




dj= vf-OL^Vj. 


Then it follows that 






4 

h = S c 'J a J ■ 
y'=i 



(68) 
(69) 

(70) 

The form of the Hermitian, 4x4 matrix, (c,y), is simplified somewhat 
by orthogonality relations: 

Cl2 = C 2 1 = C34 = C43 = 0. (71) 

Thus, by solving (70), we obtain the four coupling coefficients, a,. 

VII. RAY TRACING 

We envision that our results will find general applicability in prob- 
lems where "ray-tracing" is to be performed through anisotropic/ 
gyrotropic media and where the relative attenuation of alternative 
paths is to be explored. We shall indicate how to use the preceding 
results to this end. 

As indicated in the discussion following equation (50), there are 
difficulties in working directly with the Poynting vector, S. Even so it 
is the more macroscopic variable as it is the ray direction, as opposed 
to k which depends upon the wavefront normal. Most of the difficulty 
is resolved if we make the assumption (which does not seem too 
restrictive in practice) that the initial medium of the ray is to be 
isotropic (such as air, glass, index matching fluid, etc.). 

For isotropic media, the k direction is identical with the ray direc- 
tion. It is relatively straightforward to calculate y\\ at the first interface 
when the initial medium is isotropic. By solving (62) and, subsequently, 
calculating the corresponding S, and, if desired, the amplitudes, at, we 
can decide which of the scattered ray(s) we may wish to trace further. 

After we have chosen which scattered wave to follow, we use the 
Poynting vector to trace the ray to the next interface. We know the k 
of this ray and are not faced with the awkward task of reconstructing 
it from other information. If we are interested in exact phase relation- 
ships, these are easy to obtain: 

^(second interface) = exp[i(k-Ar)]i'(first interface). (72) 

The spatial displacement from the first to the second interface has 
been designated Ar. 

We can calculate the yj at the second interface and proceed to solve 
(62). We calculate the new S; of the scattered waves, etc., and choose 
which one we will follow to the next interface. This process can be 
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iterated as often as necessary to trace the ray of interest through an 
arbitrary optical system made up of homogeneous media and inter- 
faces. 

VIII. SUMMARY AND CONCLUSIONS 

We have developed a general formalism suitable for computation in 
the analysis of optical propagation through a succession of anisotropic, 
gyrotropic, homogeneous media. Such problems arise in a large variety 
of optical devices that are made of materials such as yig, LiNb03, 
rutile, and calcite. 
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APPENDIX 

Equivalence of Ray Direction and s. 

We shall develop a first-order perturbation theory of the general 
eigenproblem. The result is directly applicable to the calculation of 
group velocity. When the permeability and permittivity are assumed 
to be frequency independent, the group velocity must be parallel to 
the ray velocity. 

Suppose we wish to perturb the eigenproblem: 

Avi = XtBvi, (73) 

where A&nd Bare Hermitian. In addition, Sis positive definite so that 
a complete set of eigensolutions with real A, is expected. 
Our perturbation assumption is that the replacement of 

A^A + eA' (74) 

requires the corresponding substitutions: 

\i^\i + €K+ ... (75) 

v, -> v, + ev; + • • ■ (76) 

Therefore, 

(A + eA'Hxi + ev'i +.-.) = (X, + eK + • • .)fi(v, + evi + ■ • •). (77) 

Expanding and collecting terms with common powers of e, we obtain 
for €° the original (73); we obtain for e 1 : 

A vl + A ' v, = KB v< + KB v, . (78) 
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Taking the dot product with vf and cancelling terms: 

vM'v.— A'vf.fiv.. (79) 

Therefore, we can calculate A,' from 

A;=^i. (80) 

The group velocity is calculated by differentiating co(k) with respect 
to the components of k. 10 The to dependence is determined from an 
eigenproblem of the form: 



= ufw. (81) 



Thus, for example, the x-component of the group velocity, v gx is 
determined from (80) and (81): 
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v gx =da/dk x = — v* 



= (E*H Z - E*H y - H*E Z + H*E y )/2U = S x /U. (82) 

Thus, we have demonstrated: 

v* = -fUv r . (83) 
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